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PROGRAM SUMMARY 



Title of program: Optimal Jet Finder (OJF_014) 
Catalogue identifier: (supplied by the Publisher) 
Distribution format: (supplied by the Program Library) 
Computer: any computer with the FORTRAN 77 compiler 

Tested with: g77/Linux on Intel, Alpha and Sparc; Sun f77/Solaris (thwgs.cern.ch); 
xlf/AIX (rsplus.cern.ch); MS Fortran PowerStation 4.0/Win98 

Programming language used: FORTRAN 77 

Memory required: ~1 MB (or more, depending on the settings) 

Number of bytes in distributed program, including examples and test data: 1.39 MB 

Keywords: hadronic jets; jet finding algorithms 

Nature of physical problem 

Analysis of hadronic final states in high energy particle collision experiments often 
involves identification of hadronic jets. A large number of hadrons observed in the 
detector is reduced to a few jets by means of a jet finding algorithm. The jets are 
used in further analysis which would be difficult or impossible when applied directly 
to the hadrons. Reference [1] provides a brief introduction to the subject of jet find- 
ing algorithms and a general review of the physics of jets can be found in [2]. 

Method of solution 

The software we provide is an implementation of the so-called optimal jet definition 
(OJD). The theory of OJD was developed in [3], [4], [5]. The desired jet configu- 
ration is obtained as the one that minimizes ^.r, a certain function of the input 
particles and jet configuration. 

Restrictions on the complexity of the program 

The size of the largest data structure the program uses is (maximal number of 
particles in the input) x (maximal number of jets in the output) x 8 bytes. (For 
the standard settings <1 MB). Therefore, there is no memory restriction for any 
conceivable application for which the program was designed. 

Typical running time 

The running time depends strongly on the physical process being analyzed and the 
parameters used. For the benchmark process we studied, e"'"e~ — > W"'"W~ — > 4 jets, 
with the average number of ~ 80 particles in the input, the running time was 
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< 10~^ s on a modest PC per event with ntrics=l- (We took ntrics=10 for all events 
in the benchmark process, however for the majority of events ntries ~ 3 would suf- 
fice. For small values of ntries the global minimum of O/j may be missed in some 
fraction of events resulting in the deterioration of the precision of measurements 
based on the jet algorithm.) For a fixed number of jets the complexity of the algo- 
rithm grows linearly with the number of particles (cells) in the input, in contrast 
with other known jet finding algorithms for which this dependence is cubic. The 
reader is referred to [1] for a more detailed discussion of this issue. 

References 

[1] D. Yu. Grigoriev, E. Jankowski, F. V. Tkachov, e-print hep-ph/0301185[ to be 
published. 

[2] R. Barlow, Rep. Prog. Phys. 36, 1067 (1993). 

[3] F. V. Tkachov, Phys. Rev. Lett. 73, 2405 (1994); Erratum, 74, 2618 (1995). 
[4] F. V. Tkachov, Int. J. Mod. Phys. A12, 5411 (1997). 
[5] F. V. Tkachov, Int. J. Mod. Phys. A17, 2783 (2002). 



3 



Contents 

1 IntroductiorJ 



2 Algorithml 



2.1 Criterion for the optimal iet configuration! 



2.2 Algorithm for minimizing Q 



2.3 Formulas! 



3 Code and usage! 



3.1 Code and data structure! 



3.2 Normalization of energy unit d 



3.3 Error messaged 



3.4 Key minimization subroutine □■minimize 



3.5 User callable subroutined 



3.6 Compilation! 



3.7 Exampld 



4 Definitions of constants: oif_par.f'h! 



5 Common block definitions: oif_com.fh! 



5.1 Input of the eyent! 



5.2 Configuration of iets (output) 



6 Conclusion! 



Referenced 



1 Introduction 



This paper represents an official public release of the Optimal Jet Finder 
(OJF) library — a FORTRAN 77 implementation of the so-called optimal jet 
definition (OJD) for identification of jets in hadronic final states of particle 
collisions. For a brief introduction into the subject the reader is referred to 
[1], and a general review of the physics of jets can be found in [2]. The theory 
of OJD was developed in [3] , [4] , [5] , to which papers the reader is referred for 
a detailed discussion of the physical motivations behind and the derivation of 
OJD. The principal points of the theory are summarized below. 

(i) Calorimetric measurements with hadronic final states P must rely on ob- 
servables / (P) that possess a special "calorimetric", or C-continuity which is 
a non-perturbative generalization of the familiar IR safety (see [4] for details) 
and which guarantees a stability of / (P) against distortions of P such as 
caused by detectors. Ref. [4] pointed out C-continuous analogues for a variety 
of observables usually studied via intermediacy of jet algorithms. The funda- 
mental role of such observables is highlighted by two facts: (1) An observable 
inspired by [4] played an important role in the selection of top quark events in 
the fully hadronic channel at DO [8], [9]. (2) The Jet Energy Flow project [10] 
provides numerical evidence that C-continuous observables may indeed help 
to go beyond the intrinsic limitations of conventional procedure based on jet 
algorithms in the quest for the 1% precision level in the physics of jets. 

(ii) The proposition that the observed event P inherits information (as mea- 
sured by calorimetric detectors) from the underlying quark-and-gluon event q 
is expressed as 

/ (q) ~ / (P) for any C-continuous /. (1) 

(iii) For each parameter M on which the probability distribution tim (P) 
of the observed events P may depend, there exists an optimal observable 
/opt (P) = dM^^TTM (P) for the best possible measurement of M [7]. This is 
a reinterpretation of the Rao-Cramer inequality and the maximal likelihood 
method of mathematical statistics in terms of the method of moments. In the 
context of multi-hadron final states as "seen" by calorimetric detectors, such 
an observable is automatically C-continuous. 

(iv) If the dynamics of hadronization is such that eq. (1) holds, then good 
approximations for /opt could exist among functions that depend only on Q 
which is a parameterization of P in terms of a few pseudo-particles (jets), 
found from a condition modeled after eq. (1): 

/ (Q) ~ / (P) for ^-ny C-continuous /. (2) 



5 



This simply translates the meaning of jet finding as an inversion of hadroniza- 
tion into the language of C-continuous observables. 

(v) C-continuous observables can be approximated by sums of products of 
simplest such observables that are linear in particles' energies 



where a runs over all particles in the event and pa (a unit vector) denotes the 
direction of the 3-momentum of the a-th particle; / is any continuous function 
of a direction only. (The relevant theorems can be found in refs. [4] and [5].) 

(vi) So it is sufficient to explore the criterion (2) with only /'s of the form (3). 
Then one can perform a Taylor expansion in angular variables and obtain a 
factorized bound of the form 



where all the dependence on / is localized within C/^_r (so the bound remains 
vaUd for any C-contimious /) and [P, Q] is a function of the jet configu- 
ration Q (and the event P) only (the explicit form is given in section 2.1). 

(vii) Since the collection of values of all / on a given event P is naturally inter- 
preted as the event's physical information content, the bound (4) means that 
the distortion of such content in the transition from P to Q can be controlled 
by a single function; so the loss of physical information in the transition is min- 
imized if Q corresponds to the global minimum of Qr. The Optimal Jet Defi- 
nition amounts to finding the jet configuration Q which minimizes Qr [P, Q], 
depending on specific application, either with a given number of jets or with 
a minimum number of jets while satisfying the restriction fl^ [P, Q] < cUcut 
with some parameter a;cut > which is similar to the jet resolution i/cut of 
recombination algorithms. 

The beta versions of the OJF library have been available since 1999 (see 
[6]). The publicly available code [11] was first developed in the programming 
language Component Pascal [12] from the highly regarded Pascal/Modula- 
2/Oberon family, featuring a unique combination of intellectual manageability, 
safety, and efficiency of compiled code (for more on this see [13]). This made 
possible the experimentation needed to find a working and correct algorithm. 
Only after that the final port to FORTRAN 77 was performed. A subsequent 
testing on a total of O{10^^) events [14] and a substantially independent veri- 
fication [15] revealed no defects of significance, indicating a high reliability of 
the code as a result of the adopted development process. 

The OJF library can be used to obtain OJD implementations adapted for 



f{P)=EaEaf{Pa), 



(3) 



f{P)-fm<Cf,Rxnn[P,Q], 



(4) 
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specific applications (see below). 



A typical program based on OJF takes as an input the collection of particles 
(or detector cells), characterizing them by energies and directions. Then the 
optimal jet configuration is found by minimizing some function of the particles 
and jet configuration. Each particle can be shared among jets (in conventional 
schemes only an entire particle can belong to a single jet), so for an input 
event containing n particles and resulting in jets, the number of degrees of 
freedom is n- A^. A typical application (for instance, analysis of data for LHC) 
may require n ~ 100... 400 and A^ ~ 6 and therefore n - N O (1000). A large 
number of degrees of freedom makes the minimization problem non-trivial. 
The program we describe implements a practical algorithm (taking a fraction 
of a second per event on a modest personal computer) solving the optimization 
problem inherent in OJD. A key subroutine performing the minimization is 
Qjninimize. 

For a fixed number of jets the computational complexity of the algorithm 
grows as the first power of the number of particles in the input, O (n) (see 
[1] for the corresponding empirical data; for comparison, the computational 
complexity is O {'n?) for the kT algorithm [16], [17].) This feature makes OJF 
especially attractive for processing events with many calorimetric cells in the 
input. For a benchmark process we studied, e"'"e~ — > W"'"W~ 4 jets, OJF 
started to become more time efficient than the standard kT roughly at 90 cells 
in the input. ^ 

The version of OJF being described is called OJF_014. The FORTRAN 77 
source code is available from http:/ /www.inr.ac.ru/~ ftkachov/projects/jets/ 

OJF contains many interface subroutines (usually with names starting with 
set_ or get_) allowing to input data, read output and change the algorithm 
parameters. The user is not supposed to write data directly to common blocks 
but is to use the interface subroutines. Similarly, only the parameters that can 
be accessed by the user callable subroutines are supposed to be changed (at 
least at level of standard application). 

OJF has been compiled on various platforms: g77/Linux on Intel, Alpha and 
Sparc, Sun f77/Solaris (thwgs.cern.ch), xlf/AIX (rsplus.cern.ch), MS Fortran 
PowerStation 4.0/Win98. Its numerical stability was tested on a total of ~ 10^*^ 
events. 

OJF contains options to handle both the center of mass (spherical) kinematics 
of lepton collisions and the cylindrical kinematics of hadron collisions. 



^ For nntries = Ij See section 2.1. 
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2 Algorithm 



2. 1 Criterion for the optimal jet configuration 



The input for the program is an event: a coUection of n particles from the 
detector (or n hit detector cells), indexed with a = 1, 2, 3, n. Each particle 
is characterized by its energy, and its direction described by the standard 
angles 9 a, or equivalently by transverse energy, E"^, pseudorapidity, 77a, and 
the angle ipa- The a-th particle in the input is assigned the 4-momentum p^. 

Pa^ Ea- {I, sin 9a COS </?„, siu 9a sin (pa, COS 9a) (5) 



or 

Pa = ■ (cosh 77a, COS (fa, siu ipa, siuh Tja) ■ (6) 

depending on what parameters describe the particles. The output of the pro- 
gram is a set of N jets, indexed with j = 1, 2, 3, iV. The jet configuration 
is described by recombination matrix {zaj} components of which satisfy: 

< Zaj < 1 for all a,j, (7) 

N 

Zaj < 1 for all a. (8) 

The number Zaj gives the fraction of the a-th particle which goes into forma- 
tion of the j'-th jet. Each Zaj can take any value between and 1. The final 
value of the recombination matrix {zaj} is the result of the algorithm. The 
4-momentum qj of the j-th jet is defined as 

n 

Qj = Yl ^ajPa- (9) 
o=l 

The final (optimal) jet configuration is the one that minimizes the value of 
some function Q {{zaj}) depending on the recombination matrix {Zaj} and all 
Pa as parameters. 

The definition of fl follows some intermediate sub-definitions. Part of the a-th 
particle that do not go into formation of any jet: 

N 

Za = l-Y.Zaj. (10) 
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The rest of the definitions are given separately for spherical kinematics (lepton 
collisions) and for cylindrical kinematics (hadron collisions). 



Spherical kinematics. Overall energy left outside jets Esoit, called soft 
energy: 

n 

Esoh = '^ZaEa- (11) 

a=l 



The function Y, called fuzziness: 

N 

Y^2Y,q,qj, (12) 



where qj is light-like {q/ = 0) 4-direction defined: 

qj = (1, sin^jC0S</7j, sin 9 j simp j, cosOj) , (13) 



with 



cosOj = ^=M^===, (14) 
cos^,^^=^^k==, (15) 



Cylindrical kinematics. The soft energy is the overall transverse energy 
left outside the jets 

n 

E,oft = Y.^aEi. (17) 

a=l 



The fuzziness Y is defined again by (12) with qj, light-hke {qj^ — 0) 4-direction 
given by: 

qj = (cosh r]j , cos ipj , sin ipj , sinh j ) , (18) 
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where 



YZ=\ ^ajEg rja 



(19) 



and cosipj, simpj given by (15), (16). 

Finally, in both cases, ^2 is a linear combination of Y and £^soft with the 
parameter R weighting their relative contribution: 

n{{zaj}) = ^Y + E,,f,. (20) 



If we intend to reconstruct an event to some fixed number of jets the procedure 
is as follows. Start with some initial value of the recombination matrix Zaj, for 
example chosen randomly and minimize the function Q, {{zaj}) with respect 
to {zaj}- The algorithm described in section 2.2 allows one to find a local 
minimum while starting with some initial value of {zaj}- Repeat the procedure 
a few times (parameter ntries); starting each time with different initial value 
of {zaj} and take the smallest of the minima obtained at each try. The value 
of the recombination matrix {zaj} that gives the smallest of the minima is the 
final jet configuration. If the initial value of {zaj} is not chosen randomly it is 
useless to do the minimization procedure more than once as the minimization 
algorithm is deterministic. 

If the number of jets is to be determined in the process of jet reconstruction 
we can repeat the procedure described above for different number of jets 
each time. This means, we start with some minimal number of jets, e.g. = 1 
and find the corresponding A^-jet configuration and check whether 

^ iUaj}) < OOcnt (21) 



is fulfilled for that configuration. If so, this is the final jet configuration we look 
for. If not, we increase by one and repeat the A^-jet procedure for the new 
A^. We continue until (21) is finally satisfied (which has to be true for some 
sufficiently large number of jets). The parameter a;cut is some (small) positive 
number one chooses and it is analogous to the jet resolution parameter of 
conventional recombination algorithms. 



2.2 Algorithm for minimizing Q 



The domain of the function fl {{zaj}) is a (n • A^)-dimensional product of sim- 
plices. That is, for fixed a the numbers Zaj, j = 1, N, satisfying conditions 
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(7), (8) define an A^-dimensional simplex. In typical application n ~ 200 (or 
more) and ~ 5 and therefore f2 is a function of ~ 1000 variables. The 
algorithm described below allows for efficient minimization of Q {{zaj})- 

The algorithm iteratively descends into local minimum of ^}{{zaj}) starting 
from a given initial value of {zaj}- At each iteration, subsequently for each 
particle, {zaj} is moved into new position in the domain that gives the smaller 
value of fl. The iteration loop is terminated when no particle is moved at a 
single iteration, meaning that the local minimum has been found. (Or some 
safe number of maximal iterations has been exceeded.) 

We describe now in detail how {zaj} is moved in a single iteration step for a 
given particle. Denote z = (2:1, Z2, zn) with zj = Zaj and fl{z) = fl {{zaj}) 
with fixed a in both definitions. The change in Q, when we change z to z + rd 
can be described locally as 

n{z + Td) = n{z) +Ti -d + O (r^) , (22) 



where f = (/i, /at) , fj = dn (z) /dzj,i-d = Y.f=i fjdj and d = {di, d 



N) 



describes some direction. If z were not constrained to the simplex we could 
take d = — f and some r > to decrease fl. But choosing r and d we have to 
ensure that z + rd is within the simplex. Rewrite 



N N 



f-d = $:M = EM + M (23) 

j=l 3=1 



with the following definitions 

f, = f,-fj, fo = -fj: (24) 

N 

do = -Y.dj, (25) 



where J is any of 1, for which zj > (there always must be such J). 
Now d can be chosen as follows 



dj = i 



max (0, —fj) for all j — 0, A", for which Zj — 

(26) 

—fj for all j — 0, j ^ J for which Zj > 



and dj is chosen so that (25) is satisfied. With such choice of d and the proper 
parameter r the new candidate minimum z + rd will belong to the simplex 
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and n (z + rd) < Q(z) . In the above prescription the choice of J is arbitrary. 
We found it advantageous to choose J {zj > 0) such that the norm 

|(d,do)| =max{|d,| : j = 0, 1, TV} (27) 



is maximaL The choice of step length r is determined by the experimental 
finding that the minimum tends to be located at the boundary of the simplex. 
We find 

r = min (^|-^ ■ i = 0' > and dj < o| j (28) 



from the requirement that the new candidate minimum z + rd should be 
located at the boundary of the simplex and if this results in an increase of 
the value of fl then r is iteratively divided by a constant factor (~ 3) until 
minimum is found. 

An important technical implementation detail is the so-called "snapping". If 
some Zaj is small enough (i.e. z is close enough to a boundary of the simplex) 
then it is set to zero. A similar snapping is used for the direction d. 



2. 3 Formulas 



We give here the explicit formulas for derivatives fj = dVL{{zaj}) jdz^j used 
within the algorithm (derived from the definitions given in the previous sec- 
tions) . 



Spherical kinematics: 

= IVa^j - Ea. (29) 

Cylindrical kinematics: 

fj = 2p„g,- - ^ " {r]a - m) (fj si"^ % " ^^sh 77,) - E^. (30) 

3 Code and usage 



We describe a FORTRAN 77 implementation of the algorithm explained in 
the previous section. 
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3.1 Code and data structure 



The code (file ojf_014.f ) consists of subroutines (and functions) which can be 
divided in three logical groups: (i) interface subroutines, (ii) core subroutines 
and (iii) example jet search or result printing subroutines. In addition block 
data ojf -lock contains default values of some program parameters. The inter- 
face subroutines allow the user to enter input data, to read output or already 
entered input, to set or change program parameters and to obtain information 
about current program parameters. All parameters that are supposed to he set 
or changed by the user can he accessed hy these suhroutines. The same applies 
to all input and output data. The user is not supposed to write directly to com- 
mon blocks. The core subroutines (functions) perform {{zaj}) minimization 
and conversion between various data forms. The user is not supposed to call 
them directly except for Q_minimize. The subroutine Q_search is an example 
application of OJF frame to simple jet search (see section 3.5.7). The user 
may want to modify it or write their own subroutines if needed. 

All floating point variables within the program are deflned as 
DOUBLE PRECISION. If the user employs REAL type variables they should ensure 
that a proper conversion of the parameter values is made in the calls of the 
OJF subroutines. 

The file ojf_com.fh contains common block definitions of internal data struc- 
tures for OJF, for instance, matrices for parameters of input particles, output 
jets parameters and recombination matrix {zaj} ■ The file o j f _par . f h contains 
the definitions of constants used within the program. The file ojf Jkin.fh con- 
tains the definitions of two constants: sphere=l , cylinder=2. The file can be 
contained in user programs whenever reference to kinematics type is made, e.g. 



INCLUDE 'ojf_kin.fh' 
INTEGER kinematics 

kinematics=sphere 
event_setup_begin( kinematics ) 



The other two files (ojf_com.fh and ojf_kin.fh) normally do not need to be 
contained in user programs. 
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3.2 Normalization of energy units 



The energies £■„ or E-^ of input particles and the corresponding 4-momenta 
Pa are normahzed (after being entered) according to 

w ('\^ \ 

~^ Er> Pa ~^ ^ W-Lj 

Z^a=l ^a Z^a=l ^a 



for spherical kinematics or according to 

^a ~^ Pa ~^ 7=^7* FT ("^^) 

2^a=l -^a 2^a=l ^a 



for cylindrical kinematics. The normahzation constant Z)a=i -^a or J2a=i 
is stored to interpret properly the final output. The normalization allows to 
make the implementation independent of energy units and scale. 



3.3 Error messages 



Significant part of the code consists of various checks. For example: 

IF (.NOT. ojf_event_begin) THEN 

WRITE (6,*) ' add_particle : 20: wrong call sequence' 

WRITE (6,*) 'call event_setup_begin first' 

STOP 'add.particle: 20' 
END IF 

The checks are used to assure that subroutines are not called in inappropriate 
order, chosen parameters or input data do not have pathological values and 
that the program runs properly. The check can generate an error message 
and terminate the program. Messages with numbers 20-29 are due to the user 
errors. Messages with numbers > 30 are generated by program failures, so 
should you get such a message, please inform the authors; please include the 
corresponding event in text form. 



3.4 Key minimization subroutine Q-minimize 



The subroutine Q_minimize minimizes Q{{zaj}) for a given number of jets 
starting from the existing configuration of {zaj} ■ An example program that 
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uses Qjninimize is given in section 3.7. 

The subroutine performs iteratively the minimization algorithm described in 
section 2.2. The iteration loop is terminated when no particle is moved in a 
single iteration or the maximal number of iterations is exceeded. (We regard 
that the minimum is found only in the former case.) Default value of the 
maximal number of iterations is set 1000 which corresponds to ~ 1 second of 
computing time on a modest computer. It can be changed with setjnaxiter 
( mcLxiter ) , see section 3.5.3. In each iteration, a loop over all particles is 
run (a = 1, 2, n). For each particle separately, new candidate {zaj} for the 
minimum is found. The direction, d, and step , r, are computed according to 
the procedure described in section 2.2. Unless the step is zero or "infinity", 
indicating that the particle should not be moved, the condition 

n(z + Td)<n (z) . (33) 

is checked. If the condition is met the recombination matrix {zaj} is moved 
into the new position. If not, the step is reduced 3 times and (33) is checked 
again. If it is not true r is reduced again and so on. If r falls below some small 
parameter (r |(d, do)\ < eps_dist) the particle is not moved and the program 
proceeds to the next particle. 

3.5 User callable subroutines 

We describe all user callable subroutines other than Qjninimize explained 
above. 

3.5.1 Event setup 

event_setup_begin ( kinematics ) 

input: 

INTEGER kinematics kinematics type 
The subroutine begins initialization of a new event. It must be called before 
event data is entered. The parameter kinematics informs the program what 
type of kinematics is used: spherical (center of mass collisions), kinematics=l 
or cylindrical (hadron collisions), kinematics=2. If the file ojf_kin.fh is 
included constants sphere and cylinder can be used to assign value to 
kinematics: 

INCLUDE 'ojf_kin.fh' 
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INTEGER kinematics 

kinematics=sphere 
event_setup_begin( kinematics ) 

Kinematics ought to be set once for all events in a job. 

add_particle ( energy, theta, phi ) 

input: 

DOUBLE PRECISION energy energy Ea 

DOUBLE PRECISION theta angle Oa 

DOUBLE PRECISION phi angle 0„ 
The subroutine is used to enter input data. It must be called between 
event_setup_begin and event _setup_end. Each call adds a particle (=detec- 
tor cell) to the event. The energy Ea of the particle can be in any units, for 
example GeV. The direction of the particle is described by the standard angles 
9a (measured from beam axis) and (pa- 

add_particle_raw ( px, py, pz ) 

input: 

DOUBLE PRECISION px, py, pz 3-momentum components 
The subroutine is used to enter input data, as an alternative to add_particle. 
It must be called between event_setup_begin and event_setup_end. Each 
call adds a particle (^detector cell) to the event. The parameters px, py, pz 
are 3-momentum components in the same units as energy in add_particle. 
The beam axis is in z-direction. The subroutine is useful with output of Monte 
Carlo event generators. It can be freely mixed with add_particle. 

event_setup_end 

must be called after all input particles are entered and before the jet search 
can be undertaken. No particles can be added to the event afterwards. This 
subroutine is needed for internal housekeeping. For instance, it provides the 
proper normalization of the energies of the particles. 
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3.5.2 Setup of initial jet configuration 

jets_setup_begirL ( njets, Radius ) 

input: 

INTEGER njets number of jets, 

DOUBLE PRECISION Radius parameter R 
The subroutine has to be called to begin setup of the initial jet configuration - 
the initial value of the recombination matrix {zaj}, necessary for the iterative 
minimization of Q, as explained in section 2.2. It is called automatically by 
Q_search but it must be called explicitly if Q_minimize is used instead. The 
number of jets, N, must be positive. (The event will be reconstructed to the 
number of jets entered here.) R is the parameter in equation (20). It has to 
be positive and not too close to zero. The bigger R is, the less energy is left 
outside the jets. New configurations of jets can be set up any number of times 
for the same event. The value of the seed from which the random number 
generator will start for this jet configuration is stored at this point. (Prom 
this point until the first invocation of anything random, the seed can be reset 
by set_seed.) If you need only to change R and proceed with minimization 
starting from the current configuration, use reset_Radius. 

set_seed ( seed ) 

input: 

INTEGER seed 

This is to allow variation in random initial configurations of jets in case 
there are several local minima. It may be called once for a whole sequence 
of events - each event starts with a seed set up by the internal random 
number generator. The seed can be read (see get_seed) and used as a key 
to regenerate the corresponding configuration of jets (i.e. local minimum; so 
the local minimum is completely determined by its seed). It must be called 
after jets_setup_begin but cannot be called after the first invocation of 
init_z_random or init_random_all or jets_setup_end and until the next 
j ets_setup_begin. 

reset_Radius ( Radius ) 

input: 

DOUBLE PRECISION Radius parameter R 
The subroutine changes the value of the parameter R in equation (20). R has 
to be positive and not too close to zero. A large value of R means less energy 
is left outside jets. The subroutine can be called at any time - the current con- 
figuration of jets is not affected (only O, is recalculated properly). This may 
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be useful for setting up interesting variations of the algorithm ("annealing") 
in which one starts from some small value of R and then changes it gradually, 
fine-tuning the resulting jet configurations by calls to Qjninimize. With in- 
finitesimal values of R, the global minimum occurs for jet configurations with 
the most energetic particles playing the role of jets, so this can be used to 
obtain the most energetic (narrow clusters of) particles. 

init_z_random_all 

The subroutine can only be called between jets_setup_begin and 
jets_setup_end. It is the simplest way to initialize the recombination ma- 
trix {zaj}: completely and uniformly random in the direct product of all the 
simpliccs corresponding to particles. If only specific particles need to be ran- 
domized, irLit_z_random( a ) should be used. If only specific particles need 
to be set non-randomly, init_z ( a, z_in) or assign_to_jet ( a, j ) should 
be called for those particles. Then init_z_random_all can be called to ran- 
domize the remaining particles. If this is not called explicitly, the particles 
not explicitly initialized are set to "neutral" positions (democratically shared 
between all jets and the soft energy). 

assign_to_jet ( a, j ) 

input: 

INTEGER a index of the particle 

INTEGER j index of the jet 
The subroutine can only be called between jets_setup_begin and 
jets_setup_end. It can be used to set the initial configuration of jets explicitly, 

for instance, when the output of another jet algorithm is to be fine-tuned. It 
sets the value Zaj = 1 for the given a, j, i.e. directly assigns the a-th particle 
to the j-th jet. It must have 1 < a < n and < j < A^; j = corresponds 
to soft energy. The subroutine only sets the initial configuration. No elements 
of the recombination matrix are protected from being changed by subsequent 
minimizations. 

init_z_f rom ( a, z_in ) 

input: 

INTEGER a index of the particle 

DOUBLE PRECISION z_in(0 :njets_inax) components Zao, Zai, ZaN 
The subroutine can only be called between jets_setup_begin and 
jets_setup_end. It can be used to set the initial configuration of jets explicitly, 
for instance, to fine-tune the output of another jet algorithm. It initializes the 



18 



recombination matrix {zaj} for the a-th particle, i.e. sets Zao, Zai, ZaN- Only 
A^+1 components of the vector z_in(0 :njets_max) are used. The components 
must be all non-negative but do not need to be normalized correctly - correct 
normalization will be imposed automatically; z_in(0) is the particle's fraction 
relegated to soft energy. For instance, z_in( j ) can be some measure of distance 
between the a-th particle and the j'-th jet from another jet algorithm. The 
subroutine only sets the initial configuration. No elements of the recombination 
matrix are protected from being changed by subsequent minimizations. 

init_z_randoin ( a ) 

input: 

INTEGER a index of the particle 
The subroutine can only be called between jets_setup_begin and 
jets_setup_end. It does random initialization of the recombination matrix 
{zaj} for the a-th particle. 

jets_setup_end 

The subroutine must be called prior to minimization. It does housekeeping 
such as initialization of the particles whose recombination matrix elements 
have not been explicitly initialized by calls from init_z_random. 

3.5.3 Setting algorithm control parameters 

set-maxiter ( maxiter ) 

input: 

INTEGER maxiter maximal number of iterations 
The subroutine can be called to change the maximal number of iteration, see 
3.4. Default value of the maximal number of iterations is set to 1000 which 
corresponds to ~ 1 second of computing time on a modest PC. It can be called 
at any time. 

set_njets_limits ( nstart, nstop ) 

input: 

INTEGER nstart starting number of jets 

INTEGER nstop maximal number of jets 
The subroutine is needed in conjunction with Q_search only. It sets the start- 
ing and the final number of jets in Q_search (see the end of section 2.2 
and description of Q_search in section 3.5.7). The parameters must obey 
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1 < nstart < nstop and nstop < njets_mcLX (constant njetsjnax, set in 
ojf_par.fh, defines the dimension of matrices and is the maximal allowed 
number of jets). The default values are: nstart=l and nstop= njets_max=20. 
The subroutine can be called any time. 



set_ntries ( n ) 

input: 

INTEGER n number of tries 
The subroutine is needed in conjunction with Q_search only. It sets the num- 
ber of tries to find the minimum with different random initial configurations 
for each number of jets (see the end of section 2.2 and description of Q_search, 
section 3.5.7). The parameter n must be positive. The larger n, the higher the 
probability that the found configuration is the global minimum. Note that 
number of local minima correlates positively with number of hard partons. 
Usually values ~ 10 should suffice. The subroutine can be called at any time. 



set_trace_runoved ( bool ) 

input: 

LOGICAL bool see text 
The subroutine with parameter bool=.TRUE. turns on the option in which 
Q_minimize prints how many particles were moved at each iteration; with 
bool=. FALSE, it switches the option off (default). The subroutine can be 
called any time. 



3.5.4 -Access to parameters 

get-kinematics ( kinematics ) 

output: 

INTEGER kinematics type of kinematics 
The subroutine returns the type of kinematics. The possible values are 1 
(spherical kinematics, center of mass collisions) and 2 (cylindrical kinematics, 
hadron collisions), which is equivalent to constants sphere and cylinder if 
the header file ojf_kin.fh is included (see also section 3.5.1). The subroutine 
cannot be called prior to the very first call of event_setup_begin. 
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get_nparts ( nparts, e_scale ) 



output: 

INTEGER nparts number of particles in the event 

DOUBLE PRECISION e_scale total energy of the event 
The subroutine returns the number of particles, n, and the total energy in 
the event. The total energy is the sum of the usual energies of the particles 
for spherical kinematics Zla=i and the sum of transverse energies for cylin- 
drical kinematics J2a=i physical units, i.e. prior to the normalization 
J2a=iEa = 1 or J2a=iEa — ^- Other words, e_scale is the normalization 
constant. Energy/momentum parameters returned by some other subroutines 
are normalized by the value of e_scale. The subroutine cannot be called be- 
tween event_setup_begin and event_setup_end. 



get_particle ( a, e, xta, phi, p, ephys, pphys ) 



input: 
INTEGER 



index of the particle 



output: 

DOUBLE PRECISION e 

DOUBLE PRECISION xta 

DOUBLE PRECISION phi 

DOUBLE PRECISION p(0:3) 



normalized energy Ea or 
angle 9a or pseudorapidity ria 

angle 0^ 

normalized 4-momentum Pa 



DOUBLE PRECISION ephys energy E^ or E^ not normalized 

DOUBLE PRECISION pphys (0:3) 4-momentum pa not normalized 
The subroutine returns parameters of the a-th particle. For spherical kine- 
matics the parameters are the usual energy, Ea, and the standard angles 6a 
(from the beam axis) and (pa- For cylindrical kinematics the parameters are 
the transverse energy, E^, pseudorapidity, rja, and the angle (pa- The value 
of e is normalized and the ephys is in the same units as used in input, i.e. 
ephys = e ■ e_scale (see the previous subroutine for e_scale). All angles are 
in degrees. In both kinematics, p(0:3) and pphys (0:3) are the normalized 
and non- normalized 4-momenta, Pa, of the particle. The subroutine cannot be 
called between event_setup_begin and event_setup_end. 



get_njets ( njets ) 

output: 

INTEGER njets number of jets 



21 



The subroutine returns the number of jets in the current configuration of jets. 
It cannot be called before the first configuration of jets is setup. 

get_seed ( seed ) 

input: 

INTEGER seed seed for random generator 
The subroutine returns the value of the seed for the random generator, used 
for setting up the current random jet configuration. The value of the seed is 
"locked" (causing attempts to reset it to result in program termination) by 
the first invocation of anything "random" and retained until "unlocked" and 
reset by jets_setup_begin. 

get_Radius ( R ) 
output: 

DOUBLE PRECISION R parameter R in eq. (20) 
The subroutine returns current value of the parameter R in equation (20). 

getjncLxiter ( meixiter ) 

output: 

INTEGER maxiter maximal number of iteration 
The subroutine returns the maximal number of iterations (see section 3.4). 

get_njets_limits ( nstart, nstop ) 

output: 

INTEGER nstart starting number of jets 

INTEGER nstop maximal number of jets 

The subroutine returns the current values of the starting number of jets and 
the maximal number of jets in the subroutine Q_search (see the end of section 
2.2 and description of Q_search in section 3.5.7). 

get_ntries ( n ) 

output: 

INTEGER n number of tries 
The subroutine returns the current number of tries in Q_search, the number 
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of attempts to find minimum with different random initial configurations for 
each number of jets (see the end of section 2.2 and description of Q_search, 
section 3.5.7). 

3.5.5 Access to results 

get_criterion ( omega, y, esoft ) 

output: 

DOUBLE PRECISION omega value of Vt 

DOUBLE PRECISION y value of Y 

DOUBLE PRECISION esoft value of E^oit 
The subroutine returns the value of y and £^soft- Whenever a jet configu- 
ration is set up or modified, the corresponding values of Vt, Y and £^soft ^-re 
recalculated and can be retrieved using this subroutine. 

get_jet ( j, e, xta, phi, q, qtilde, ephys, qphys ) 
input: 

INTEGER j index of the jet 

output: 



DOUBLE 


PRECISION 


e 


normalized energy 








or normalized transverse energy 


DOUBLE 


PRECISION 


xta 


angle 6j or pseudorapidity rjj 


DOUBLE 


PRECISION 


phi 


angle (pj 


DOUBLE 


PRECISION 


q(0:3) 


normalized 4-momentum qj 


DOUBLE 


PRECISION 


qtilde (0:3) 


4-direction qj 


DOUBLE 


PRECISION 


ephys 


energy (or transverse energy) 








in physical units 


DOUBLE 


PRECISION 


qphys (0:3) 


4-momentum in physical units 



The subroutine returns parameters of the j-th jet, where j obeys < j < N 
and j = is the zeroth "jet" , name for the fractions of particles that do not 
belong to any jet (i.e. soft energy). For spherical kinematics the parameters 
are the usual energy, Ej, normalized e and non-normalized ephys (i.e. in the 
units of energy used in the input), the standard angles 6j (from the beam axis) 
and (pj. For cylindrical kinematics the parameters are the transverse energy, 
E^, normalized e and non-normalized ephys (i.e. in the units of energy used 
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in the input), pseudorapidity rjj and the standard angle (f)j. All angles are in 
degrees. For both kinematics the parameters q(0 : 3) and qtilde (0 : 3) are the 
normalized and non-normalized 4-momentum of the jet. qj is the 4-direction 
defined in section 2.1. 



get_z ( a, z_out ) 



input: 




INTEGER a 


index of the particle 


output: 




DOUBLE PRECISION z_out (0 :njets_max) 


components ZaO,Zal,..., ZaN 



The subroutine returns the components Zao, Zai, ZaN of the recombination 
matrix for the a-th particle, a must satisfy 1 < a < n. The value of z_out ( j ) 
is the a-th particle contribution to the j-jet and j = corresponds to the soft 
energy. Note: E5=o Z-Oiit(j) = 1- 



get_particle_split ( a, total_jets, jet, zj ) 

input: 

INTEGER a index of the particle 

output: 

INTEGER total_jets number of jets the particle 

belongs to 

INTEGER j et (0 : nj ets_max) indices of the jets 

DOUBLE PRECISION zj (0 :njets_max) corresponding ^r^j 
For the a-th particle the subroutine returns: the number of jets (including soft 
energy 0-th "jet") which include a non-zero fraction of the particle {zaj ^ 0), 
the labels of the jets and the corresponding values of Zaj in such an order 
that zj (k) > zj (j + 1). In other words, the vector zj (0 :njets_max) is the 
collection of Zaa^Zax^ ■■■,ZaN ordered by their value (descending from the left 
to right); only the components (O:total_jets-1) are different from zero. 
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get_jet_split ( j, nwhole, whole_a, nfract, fract_a, fract_z ) 



input: 
INTEGER 



index of the jet 



output: 
INTEGER 

INTEGER 
INTEGER 

INTEGER 



nwhole number of particles entirely 

belonging to the jet 

whole_a(0 :nparts_maj?;) labels of the particles 

nfract number of particles 

belonging in some fraction 

f ract_a(0 :nparts_max) labels of the particles 



DOUBLE PRECISION f ract_z(0 :nparts_max) the fraction ^, 



aj 



For the j-th jet (0 < j < A^, j = is the soft energy) the subroutine returns: 



• number of the particles wholly in the jet, i.e. Zaj = 1 

• vector whole_a(0 :nparts_max) with labels of such particles (indices a) 

• number of particles partially in the jet, i.e < Zaj < 1 

• vector f ract_a(0 :nparts_max) with labels of such particles (indices a) 

• vector f ract_z(0 :nparts_max) with corresponding Zaj for such particles. 

The latter two vectors are synchronously ordered so that subsequent compo- 
nents of f ract_z(0 :nparts_max) do not increase. 



3.5.6 Sample print routines 
print_z_raw 

is an example subroutine to print the recombination matrix {zaj}. A possible 
output may look like: 





background 




1 




2 




3 


1 


0.0000 


0. 


,0000 


0, 


,0000 


1, 


.0000 


2 


0.0000 


0, 


,0000 


0, 


.0000 


1, 


.0000 


3 


0.0000 


0, 


,0000 


0, 


.0000 


1, 


.0000 


4 


0.0000 


0, 


,0000 


0, 


.0000 


1, 


.0000 



print _z_nice 

is an example subroutine to print the recombination matrix {zaj}. A possible 
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output may look like: 



recombination matrix z by particle label a: 

jet numbers 



a background 12 3 



1 - - - 1. 

2 - - - 1. 

3 - - - 1. 

4 - - - 1. 



print_jets 

is an example subroutine to print properties of jets. See the output of the 
example program in section 3.7.2. 



print_particles 

is an example subroutine to print properties of particles. A possible output 
may look like: 



Configuration by particle: 



EC/,) 



theta 



(soft energy is denoted as jet=0) 
phi jet [ fraction ]; ... 



0.5100 
0.4000 
0.4000 
0.2000 



6.7194 70.0000 0.0000 3 

5.2701 90.0000 0.0000 3 

5.2701 85.0000 10.0000 3 

2.6350 84.0000 -10.0000 3 



24 
25 



0.2000 
7.0000E-02 



2.6350 
0.9223 



170.0000 
90.0000 



-7.0000 
170.0000 



TOTAL: 7.5900 



100.0000 



3.5.7 Example subroutine of straightforward jet search Qsearch 

This is a simple jet search subroutine using Qjninimize as a key component. 
It is possible that the user may want to modify it, for example, when trying to 
do something with local minima. This subroutine uses only interface routines; 
it does not access internal data. 

The subroutine tries to find the configuration of jets which minimizes fl and 
ensures that < a;cut with the minimal number of jets (njets) starting from 



26 



the number of jets previously set via set_njets_start (usually the same for 
all events). For each number of jets, the search is repeated ntries times, each 
time with a different random initial value of the recombination matrix {zaj} 
and the configuration with the lowest value of Q is retained as a result. Failure 
of the search is signaled by the condition njets=0. 

Note that Q_search randomizes the initial value of {zaj}, so it is meaningless 
to use it if one wants to specify the initial configuration for {zaj}. In this case, 
the user should use Qjninimize directly. We comment that some other control 
options could be to continue attempts until a specified number of attempts 
fails to yield a better configuration or to stop the search for new minimum 
if, for example, the first three random initial configurations yielded the same 
configurations (the event has a single local minimum which is automatically 
the global one; this is quite hkely and may be useful if CPU time is an issue) . 

3. 6 Compilation 

Optimal Jet Finder consists of the following files: 

• ojf_014.f main file contains all subroutines and functions 

• ojf_com.fh contains definitions of common blocks 

• ojf_par.fh contains definitions of parameters 

• ojf _kin.fh contains definition of kinematics type parameters 

The example programs example. f, wwl60.f, wwl60a.f with input or output 
files example, in, wwl60.in, wwl60.out, wwl60a.out are added. 

To compile and run any of example programs with OJF under Linux equipped 
with g77 the user can type: 

g77 user_program. f ojf_014.f -o executable_f ile (enter) 
executable_f ile (enter) 

where user_program. f is the name of the user own program applying OJF. 
Each example program example. f, wwl60.f or wwl60a.f can be used in its 
place. The files ojf_com.fh, ojf_par.fh, and ojfJkin.fh should be available 
in the current directory (but not compiled). 

3. 7 Example 

The simplest possible example, file example . f below, should give the idea how 
Optimal Jet Finder can be used. The file example . in contains input data. 
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Each line corresponds to one particle and consists of Ea, Oa and 0a for that 
particle. The user is encouraged to study subroutine Q_search and programs 
wwl60.f, wwl60a.f providing additional, more advanced examples. 



3. 7. 1 examp le.f 

PROGRAM simplest.example 



INCLUDE 'ojf.kin.fh' 



DOUBLE PRECISION Radius 

DOUBLE PRECISION e, theta, phi 

DOUBLE PRECISION o_fin, y_fin, e_fin 

INTEGER a, seed, nparts, njets, kinematics 

LOGICAL success 



number of jets is chosen 
njets = 3 

seed for random generation of the recombination matrix 

seed = 13 
R parameter from equation (20), section 2.2 

Radius = 1.0 
choose spherical (lepton collisions) kinematics 

kinematics = sphere 

file with input data is opened 

OPENCIO, FILE='example.in' , FORM= ' formatted ' , STATUS = 'oldO 

input event setup starts 

CALL event_setup_begin ( kinematics ) 

loop over all particles in the event 
nparts = 
DO a = 1, 1999 

READ (10,*, end=1000, err=1000) e, theta, phi 

CALL add_particle ( e, theta, phi ) 

nparts = nparts + 1 
ENDDO 



1000 CLOSE(IO) 



input of the event ends 

CALL event_setup_end 
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set up random the initial value of the recombination matrix 
CALL jets_setup_begin ( njets, Radius ) 
CALL set_seed ( seed ) 
CALL init_z_raiidom_all 
CALL jets_setup_end 

minimize fl 

CALL Q_minimize ( success ) 

IF (.NOT. success) STOP 'minimum not found' 

get and print the values of Q, Y and Esoh for the final jet configuration 
CALL get_criterion ( o_fin, y_fin, e_fin ) 

WRITE(*,*) 'Omega =', o_fin 
WRITEC*,*) 'Y y_fin 
WRITE(*,*) 'E.soft =', e_fin 

prints properties of the resulting jets 
call print_jets 

END 



3.7.2 Output of the example 

Omega = . 293404849 
Y = 0.0338528071 

E.soft = 0.259552042 

SPHERE: 3 jets processed 

Configuration by jet: 

jet E E(%) theta phi 

1 1.380 18.1818 138.3848 -52.8876 

2 1.220 16.0738 124.4338 -26.6115 

3 3.020 39.7892 81.0226 0.3566 



TOTAL: 5.6200 74.0448 

Particle content by jet: 

jet label 1(3 particle (s) ) 
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E(y„) = 18.18 

theta = 138.4 

phi = -52.89 

3 particle (s) in jet as a whole: 21 22 24 

jet label 2(4 particle(s) ) : 
EC/J = 16.07 
theta = 124.4 
phi = -26.61 

4 particle (s) in jet as a whole: 17 18 19 20 

jet label 3(8 particle(s) ) : 
E(%) = 39.79 
theta = 81 . 02 
phi = . 3566 
8 particle (s) in jet as a whole: 12 3 4 
5 6 7 8 

soft energy ( 10 particle (s) ): 

10 whole particle (s) in soft energy: 9 10 11 
12 13 14 15 16 23 25 
no particles partially in soft energy 



4 Definitions of constants: ojf_par.fh 



In this section we explain the meaning of the parameters defined in the header 
file ojf_par.fh and give their default values. 

INTEGER njets_max=50 

The maximal number of jets; used for example to define the size of matrices. 
INTEGER npartsjnax=2000 

The maximal number of particles in the event; used for example to define the 
size of matrices. 

DOUBLE PRECISION zero=0 
DOUBLE PRECISION one=l 
DOUBLE PRECISION inf=10^°° 
are the numerical constants. 

DOUBLE PRECISION eps_snap=10-^ 

If Zaj < eps_snap than Zaj is set to zero, i.e. the particle is snapped to the 
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boundary of the simplex. The parameter is used in subroutines z_snap and 
z_assert. 

DOUBLE PRECISION eps_round=10-6 
DOUBLE PRECISION eps_sum=10-^ 
DOUBLE PRECISION eps_sumO=10-^ 
DOUBLE PRECISION eps_suml=10-^ 

The constants arc used to keep control of rounding errors. If some variable 
exceeds the allowing range of values more than eps_, the error message is gen- 
erated and the program is terminated. The constants are used in the subrou- 
tines d_minus_snap, z_snap, d_assert, z_assert, z_f orce_to_simplex and in 
the function pos_prod. 

DOUBLE PRECISION eps_noriii=10-^ 

The constant is used to determine whether the norm of the 3-vector (or 
transverse part of the norm in case of cylindrical dynamics) is zero. It is used 
in the subroutine j_eval_nonlinear. 

DOUBLE PRECISION eps_Et=10~6 

The constant is used to handle small values of the transverse energy of a jet. 
It is used in the subroutines grad_Y and j_eval_nonlinear. 

DOUBLE PRECISION eps_dist=10~^ 

Sec section 3.4. The constant is used to determine when to stop subsequent re- 
ductions of the step r. The constant is used in the subroutine Q_minimize_wrt. 

DOUBLE PRECISION eps_radius=10-3 

The constant sets the limit on the smallest value of i?, the parameter from 
equation (20). It is used in the subroutines jet_setup_begin and reset_Radius. 

DOUBLE PRECISION inf _step=10^° 

Sec section 3.4. "Infinite" step means that the particle should not be moved. 
The constant is used in subroutines Q_minimize_wrt, d_eval_step and 
z_move_by. 

It is imaginable that some of the parameters above may need to be changed 
but the user is advised to be careful when doing this. In particular, smaller 
values of some parameters would enhance sensitivity to rounding errors, caus- 
ing the safety checks to generate error messages and terminate the program. 
One may change eps_snap to a smaller value, say 10~^, and see if the results 
would change; for a small fraction of events this may slow the finding of jets 
but help to better identify local minima. 

INTEGER random_m=259200 

The constant is used by the random number generator, the subroutine seed 
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and the function randomC). 

The constants below play a technical role and are not supposed to be changed. 
The reason for defining them is cleared in the next section. 

INTEGER par_Et=4 
INTEGER par_eta=5 
INTEGER par_Eteta=6 
INTEGER par_y=7 
INTEGER par_pOshmpzch=8 
INTEGER par_tilde=9 



5 Common block definitions: ojf_com.fh 

The header file ojf_com.fh contains common block definitions and is included 
in most of the subroutines. The user is not supposed to write to common blocks 
directly but to use interface subroutines. Data that cannot be accessed that 
way is not supposed to be used by the user. 



5. 1 Input of the event 



LOGICAL ojf _event_begin, ojf _event_set 

The two logical values bracket the event setup: 

FALSE, FALSE - at start of program, no event has been set up; 

TRUE , FALSE - event setup in progress, adding particles; 

FALSE, TRUE - event setup completed, can search for jets. 

INTEGER ojf _kinematics 

The variable marks the type of kinematics: 1 - spherical kinematics of lepton 
collisions, 2 - cylindrical kinematics of hadron collisions. 

INTEGER ojf_nparts 

The number of particles in the event. 



k 
k 
k 



COMMON /ojf .event/ 

ojf _event_begin, ojf _event_set , 
ojf .kinematics, ojf_nparts, 
ojf_p, ojf_e, ojf_e_scale 
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DOUBLE PRECISION ojf_p(0:6, 1 :nparts_max) 

The matrix stores the properties of particles: 



ojf_p(0, 


particle_label) 


energy K 


ojf_pCl, 


particle_label) 


x-component oi momentum 


ojf_p(2, 


particle_label) 


y-component or momentum pa 


ojf _p(3, 


particle_label) 


z-componcnt of momentum 


ojf_p(4, 


particle_label) 


transverse energy E-^ 


ojf_p(5, 


particle_label) 


pseudorapidity r]a 


ojf_p(6, 


particle_label) 


c(niil)iiiati()n E-^ ■ //„ 



and particle_label is the index a of the particle. The constants par_Et=4, 
par_eta=5, par_Eteta=6 arc defined to access the components of the matrix, 
e.g. ojf _p(par_Eteta, particle_label) . 

DOUBLE PRECISION ojf _e (1 :nparts_max) 

The vector stores the energies of the particles. 

DOUBLE PRECISION ojf_e_scale 

The variable stores the energy scaling factor (see section 3.2). 
5.2 Configuration of jets (output) 

COMMON /ojf.jets/ 
& ojf _jets_begin, ojf _jets_set , 

& ojf_njets, ojf_seed, ojf_Radius, 

& ojf_z, ojf_b, ojf_q, 

& ojf_Omega, ojf_Y, ojf_Esoft 

LOGICAL ojf _jets_begin, ojf_jets_set 
The two logical values bracket setup of initial jet configuration: 
FALSE, FALSE - at start of program, or after event set up; 
TRUE, FALSE - jets setup in progress, change anything; 
FALSE, TRUE - jets setup complete, can do minimization. 

INTEGER ojf_njets 

The number of jets in the current configuration. 
INTEGER ojf_seed 

The seed used to generate the current (random) jet configuration. 
DOUBLE PRECISION ojf_Radius 
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The value of R, the parameter in equation (20). 



DOUBLE PRECISION ojf_z(0:njets_max, 1 :nparts_max) 
The recombination matrix, {zaj}- 

LOGICAL ojf _b(0:njets_max, 1 :nparts_max) 

It is used to indicate that the particle belongs to (TRUE) or does not belong 
(FALSE) to the boundaries of the simplex, i.e. Zaj = 0. 



DOUBLE PRECISION ojf_q(0:12, Irjetsjnax) 



The matrix s 



ores the properties of particles: 





.q( 


0, 


jet. 


.label; 


energy Ej 




.q( 


1, 


jet. 


.label; 


x-component oi momentum qj 


ojf 


.q( 


2, 


jet. 


.label) 


y-component oi momentum qj 


ojf 


.q( 


3, 


jet. 


.label) 


z-component of momentum qj 


ojf 


-q( 


4, 


jet. 


.label) 


transverse energy Ej- 


ojf 


-q( 


5, 


jet. 


.label) 


pseudorapidity rjj 


ojf 


.q( 


6, 


jet. 


.label) 


combination Ej- • rjj 


ojf 


.q( 


7, 


jet. 


.label) 


fuzziness Y 


ojf 


-q( 


8, 


jet. 


.label) 


combination (gj)° ■ {qj)^ — (qj)^ ■ (qj)^ 


ojf 


-q( 


9, 


jet. 


.label) 


0-component of 4-direction qj 


ojf 


-q( 


10, 


jet. 


.label) 


x-component of 4-direction qj 


ojf 


.q( 


11, 


jet. 


.label) 


y-component of 4-direction qj 


ojf 


.q( 


12, 


jet. 


.label) 


z-component of 4-direction qj 



jet_label is the index j of the jet. The constants par_Et=4, par_eta=5, 
par_Eteta=6, par_y=7, par_p0shmpzch=8, par_tilde=9 are defined to ac- 
cess the components of the matrix, e.g. ojf_q(par_y, jet_label). 



DOUBLE PRECISION ojf_Omega, ojf_Y, ojf_Esoft 

The variables store the values of fl, Y and Egoit, see equation (20). 

The remaining common block /ojf_work/ contains the definitions of "work" 
variables, mainly of the types explained above. The "work" variables are used 
at the intermediate stages of computations. 
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6 Conclusion 



We implemented an algorithm for minimization of fi, the criterion for finding 
the optimal jet configuration. The user is provided with the example program 
and subroutine Q_search as the simplest possible scenarios of jet finding. It 
is not inconceivable that practical applications will require more sofisticated 
ways of using the provided library. To facilitate this, we supplied the user with 
numerous interface subroutines allowing for easy control of the program. 

We did not explore all possible ideas for optimizing of jet search as we wanted 
to leave the program generic. For example, in some cases it is possible to apply 
O JF in a few stages. First, apply it to particles with the highest energies, then 
use the result as the initial value of the recombination matrix for the next 
stage at which the rest of particles are added (or some next portion). It may 
increase the speed of finding jets without significant deterioration of the results 
[1]. However this depends on the particular physical process at hand. 
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